В этом практикуме я исследовала геном бактерии Klebsiella variicola, для которого в прошлом семестре получала протеом.
Идентификатор выбранного протеома — UP000789617
На странице протеома находим Taxon ID: 244366
Команда для скачивания архива протеома:
wget 'https://rest.uniprot.org/uniprotkb/stream?compressed=true&format=txt&query=proteome:UP000789617' -O UP000789617.swiss.gz
Код доступа сборки GenBank — GCA_905335775.2
Код доступа из RefSeq — GCF_905335775.2
Страница сборки в Genome
Скачаем геномную последовательность и таблицу локальных особенностей(в формате gff3):
datasets download genome accession GCF_905335775.2 --include genome,gff3
Распакуем скачанный архив:
unzip ncbi_dataset.zip
Найдем какой вариант генетического кода использует данный организм:
efetch -db taxonomy -id 244366 -format xml | xtract -pattern GeneticCode -element GCId
Получаем, стандартный генетический код — 11.
Найдем открытые рамки считывания с помощью программы getorf пакета EMBOSS:
getorf -sequence ncbi_dataset/data/GCF_905335775.2/GCF_905335775.2_AN2335v1_cp_genomic.fna -find 0 -minsize 150 -table 11 -outseq translated_seq.faa
Параметр -find 0 указывает на поиск рамки между стоп-кодонами. Параметр -minsize 150 указывает на минимальную длину последовательности.
Создадим белковую базу данных для blastp с названием proteome:
makeblastdb -in translated_seq.faa -dbtype 'prot' -out proteome
Убедимся с помощью команды infoseq, что транслированные последовательности не короче 50 аминокислот:
infoseq -sequence translated_seq.faa -only -length |sort -n | head
Получаем последовательности 3х прокариотических ДНК-метилтрансфераз:
echo -e "sw:P0AED9 sw:P0AEE8 sw:P23941" | tr ' ' '\n' | seqret -filter @stdin -outseq query.fasta
За счет @stdin seqret читает USA не из файла.
Проведем blastp по созданной ранее базе proteome с последовательностями ДНК-метилтрансфераз, полученными в прошлом пункте.
blastp -query query.fasta -db proteome -out blast_results.txt -outfmt 7
Выдача blastp
Лучшая находка имеет E-value = 0.00, что свидетельствует о том, что это гомологичная находка.
Это гомолог P0AED9 (Dcm, m5C-МТаза, E.coli). Некоторые ее параметры:
Название рамки: NZ_CAJOXS020000001.1_567
Координаты в геноме: 117221 - 118960
Вес находки (bit score): 733
Определим какие CDS из таблицы локальных особенностей рассматриваемого генома, находятся рядом.
Для начала отберем строки CDS той же последовательности, что и находка (NZ_CAJOXS020000001.1), а столбцы с координатами (4,5), цепью (7) и дополнительной информацией (9). Запишем их в файл CDS.tsv:
grep '^NZ_CAJOXS020000001.1' ncbi_dataset/data/GCF_905335775.2/genomic.gff | grep 'CDS' | cut -f4,5,7,9 > CDS.tsv
Теперь найдем именно соседей, добавив строку нашей находки и отсортировав, отобрав соседние строки:
echo -e '117221\t118960\t+\tNZ_CAJOXS020000001.1_567' | cat - CDS.tsv | sort -n -k1,2 | grep -C 3 'NZ_CAJOXS020000001.1_567' > neighbors.tsv
Выдача соседних CDS
Анализируя выдачу, я нашла CDS, которая достаточно сильно пересекается с найденной рамкой.
Ее координаты: 117563 - 118963, а код доступа: WP_012967658.1. Причем в аннотации указано, что это ДНК-метилтрансфераза (DNA cytosine methyltransferase)!
Попробуем определить CDS по аннотациям. Сначала найдем все белковые записи находки(по AC в нуклеотидной записи), затем фильтруем по аннотациям(используя код ДНК-метилтрансферазы), получаем AC белковых записей.
elink -db 'nuccore' -id 'NZ_CAJOXS020000001.1' -target 'protein' | efilter -query '2.1.1.37[EC/RN Number]' | efetch -format 'acc'
Однако программа ничего не нашла. Возможно это связано с отсутвием соотвествующих аннотаций.